
********************************************************************************
**      REPLICATION FILES FOR "PRECAUTION AND PROPORTIONALITY IN PANDEMIC     **
**      POLITICS: DEMOCRACY, STATE CAPACITY, AND COVID-19 RELATED SCHOOL      ** 
**      CLOSURES AROUND THE WORLD" BY CRONERT IN JOURNAL OF PUBLIC POLICY     **
********************************************************************************


********************************************************************************
// Setup
********************************************************************************

* Change work directory as appropriate
cd ".../JPP Replication/Dataverse"

use "JPP_covid19_2022_replication.dta", clear

* Set graph scheme
set scheme s1mono

* Set up survival analysis
stset time, failure(schoolsclosed) id(ccode) scale(1)

* Define control variables
global controls weekend ln_cases1 reg_any_case
global controls2 ln_gdpcap wdi_beds_p wdi_pop14 wdi_popurb acd_res_nuclstem
global controls3 ictd_taxinsc gd_ptss



********************************************************************************
// Summary statistics (Table S1)
********************************************************************************

sutex schoolsclosed v2x_polyarchy2019 wbgi_gee2019 fh_ipolity2 icrg_qog van_comp van_part parreg parco bureaucratic_prof bureaucratic_influence $controls $controls2 $controls3 h_polcon5 dpi_presidentialism v2x_feduni if _st & daterange, labels file("tableS1.csv") minmax replace



********************************************************************************
// Figure 1
********************************************************************************

preserve
keep if inrange(date,21942,22012)
collapse (sum) schoolsclosed any_case, by(date)
label var any_case "1+ confirmed case of COVID-19"
label var schoolsclosed "School Closures"
format date %tdMon_dd
line any_case schoolsclosed date, ytitle("Number of countries") yscale(range(0,175)) ylabel(0(25)175, grid) xscale(range(21942,22012)) xlabel(21942(14)22012) name(school_closures, replace) xtitle("") scheme(s1mono) legend(region(lstyle(none))) lcolor(black gray) lpattern(dash solid)
graph export "fig1.pdf", as(pdf) replace
restore



********************************************************************************
// Figure 2
********************************************************************************

preserve
* Omit China, Hong Kong and Mongolia
drop if country == "China"
drop if country == "Mongolia"
drop if country == "Hong Kong"

* Censor Nicaragua
replace date_firstclosure = mdy(4,7,2020) if country == "Nicaragua"

sum v2x_polyarchy2019, detail
gen median_dem = r(p50)
gen above_median_dem = 0 
recode above_median_dem (0=1) if v2x_polyarchy2019 > median_dem

gen days_since1case = date_firstclosure - date_firstcase
tab days_since1case if inlist(cabbr2,"KM")
replace days_since1case = -41 if inlist(cabbr2,"VU","LS","SB")

gen abovemed = days_since1case
label var abovemed "Above Median Democracy"
gen belowmed = days_since1case
label var belowmed "Below Median Democracy"

gen cabbr2_label = cabbr2
replace cabbr2_label = "" if days_since1case < 20 & days_since1case > -3 & wbgi_gee2019 > 1.8 & cabbr2 != "QA"
replace cabbr2_label = "" if cabbr2 == "FR"
replace cabbr2_label = "" if cabbr2 == "BE"
replace cabbr2_label = "" if cabbr2 == "SD"
replace cabbr2_label = "" if cabbr2 == "US"
replace cabbr2_label = "" if cabbr2 == "FI"
replace cabbr2_label = "" if cabbr2 == "KR"
replace cabbr2_label = "" if cabbr2 == "JP"

twoway (scatter abovemed wbgi_gee2019 if above_median_dem == 1 & date == min_date, yscale(range(-40,100)) yline(0, lcolor(black)) ytitle("Days to School Closure or" "Censoring (April 7) After First Case") mcolor(gray) msize(vsmall) mlabel(cabbr2_label) mlabcolor(gray) mlabsize(small)  legend(cols(2) size(small) region(lstyle(none)))) (qfit abovemed wbgi_gee2019 if above_median_dem == 1 & date == min_date,  yline(0) yscale(range(-40,100)) lcolor(gray)) (scatter belowmed wbgi_gee2019 if above_median_dem == 0  & date == min_date, mcolor(black) msymbol(triangle_hollow) msize(vsmall) mlabsize(small)  mlabel(cabbr2_label) mlabcolor(black)) (qfit belowmed wbgi_gee2019 if above_median_dem == 0  & date == min_date, lcolor(black) lpattern(dash))
graph export "fig2.pdf", as(pdf) replace
restore



********************************************************************************
// Models A, A+, B, B+, C, and D (Figure 3 and Table S2)
********************************************************************************

est clear
eststo A: stcox v2x_polyarchy2019 wbgi_gee2019 i.region $controls if daterange, vce(robust) strata(date_firstcase) 
di 2.71828^_b[v2x_polyarchy]
di 2.71828^_b[wbgi_gee]
estat phtest, de

eststo Ap: stcox v2x_polyarchy2019 wbgi_gee2019 i.region $controls c.v2x_polyarchy2019#c.wbgi_gee2019 if daterange, vce(robust) strata(date_firstcase)
eststo B:  stcox v2x_polyarchy2019 wbgi_gee2019 i.region $controls $controls2 if daterange, vce(robust) strata(date_firstcase)
eststo Bp: stcox v2x_polyarchy2019 wbgi_gee2019 i.region $controls $controls2 c.v2x_polyarchy2019#c.wbgi_gee2019 if daterange, vce(robust) strata(date_firstcase)

eststo C: stcox  van_comp2018 van_part2018 wbgi_gee2019 i.region $controls if daterange, vce(robust) strata(date_firstcase)
eststo D: stcox  van_comp2018 van_part2018 wbgi_gee2019 i.region $controls $controls2 if daterange, vce(robust) strata(date_firstcase)

* Table S2 (Baseline measures)
esttab A Ap B Bp C D , se(3) nogaps b(3) aic bic nonotes label star(* 0.10 ** 0.05 *** 0.01) nomtitles replace eform
esttab A Ap B Bp C D using "tableS2.csv", se(3) nogaps b(3) aic bic nonotes label star(* 0.10 ** 0.05 *** 0.01) nomtitles replace eform

* Figure 3 (Baseline measures)
coefplot (A, label(A: Baseline Model) msymbol(circle_hollow) mcolor(black) ciopts(lcolor(black))) (B, label(B: A + Additional Controls) msymbol(circle) mcolor(black) ciopts(lcolor(black))) (C, label(C: Dimensions of Democracy) msymbol(diamond_hollow) mcolor(gray) ciopts(lcolor(gray))) (D, label(D: C + Additional Controls) msymbol(diamond) mcolor(gray) ciopts(lcolor(gray))) , headings(v2x_polyarchy2019 = "{bf:Baseline Variables}" ln_gdpcap = "{bf:Additional Controls}" van_comp2018 = "{bf:Dimensions of Democracy}")  drop(*.region weekend _cons) legend(region(lstyle(none)))  xline(1) name(coefs, replace) xsize(7) ysize(4) eform(*)   scale(0.9) scheme(s1mono)
graph export "fig3.pdf", as(pdf) replace



********************************************************************************
// Models  E, E+, F, F+, G, and H (Figure S1 and Table S3)
********************************************************************************

eststo E:  stcox fh_ipolity2 icrg_qog i.region $controls if daterange, vce(robust) strata(date_firstcase)
eststo Ep: stcox fh_ipolity2 icrg_qog i.region $controls c.fh_ipolity2#c.icrg_qog if v2x_polyarchy2019 !=. & daterange, vce(robust) strata(date_firstcase)
eststo F:  stcox fh_ipolity2 icrg_qog i.region $controls $controls2 if daterange, vce(robust) strata(date_firstcase) 
eststo Fp: stcox fh_ipolity2 icrg_qog i.region $controls $controls2 c.fh_ipolity2#c.icrg_qog if v2x_polyarchy2019 !=. & daterange, vce(robust) strata(date_firstcase)

eststo G: stcox parcomp parreg icrg_qog i.region $controls if daterange, vce(robust) strata(date_firstcase)
eststo H: stcox parcomp parreg icrg_qog i.region $controls $controls2 if daterange, vce(robust) strata(date_firstcase)

* Table S3 (Alternative measures)
esttab E Ep F Fp G H , se(3) nogaps b(3) aic bic nonotes label star(* 0.10 ** 0.05 *** 0.01) nomtitles replace eform
esttab E Ep F Fp G H using "tableS3.csv", se(3) nogaps b(3) aic bic nonotes label star(* 0.10 ** 0.05 *** 0.01) nomtitles replace eform

* Figure S1 (Alternative measures)
coefplot (E, label(E: Baseline Model) msymbol(circle_hollow) mcolor(black) ciopts(lcolor(black))) (F, label(F: E + Additional Controls) msymbol(circle) mcolor(black) ciopts(lcolor(black))) (G, label(G: Dimensions of Democracy) msymbol(diamond_hollow) mcolor(gray) ciopts(lcolor(gray))) (H, label(H: G + Additional Controls) msymbol(diamond) mcolor(gray) ciopts(lcolor(gray))) , headings(fh_ipolity2 = "{bf:Baseline Variables}" ln_gdpcap = "{bf:Additional Controls}" parcomp = "{bf:Dimensions of Democracy}")  drop(*.region weekend _cons) legend(region(lstyle(none)))  xline(1) name(coefs2, replace) xsize(7) ysize(4) eform(*) scale(0.9) scheme(s1mono)
graph export "figS1.pdf", as(pdf) replace



********************************************************************************
// Models I, I+, J, J+, K, and L (Table S4)
********************************************************************************

eststo I: stcox v2x_polyarchy2019 bureaucratic_prof bureaucratic_influence i.region $controls if daterange,  vce(robust) strata(date_firstcase)
eststo Ip: stcox c.v2x_polyarchy2019##(c.bureaucratic_prof c.bureaucratic_influence) i.region $controls if daterange, vce(robust) strata(date_firstcase)
eststo J: stcox v2x_polyarchy2019 bureaucratic_prof bureaucratic_influence i.region $controls $controls2 if daterange, vce(robust) strata(date_firstcase)
eststo Jp: stcox c.v2x_polyarchy2019##(c.bureaucratic_prof c.bureaucratic_influence) i.region $controls $controls2 if daterange, vce(robust) strata(date_firstcase)

eststo K: stcox v2x_polyarchy2019 wbgi_gee2019 i.region $controls $controls3  if daterange, vce(robust) strata(date_firstcase)
eststo L: stcox v2x_polyarchy2019 wbgi_gee2019 i.region $controls $controls2 $controls3 if daterange, vce(robust) strata(date_firstcase)

esttab I Ip J Jp K L, se(3) nogaps b(3) aic bic nonotes label star(* 0.10 ** 0.05 *** 0.01) nomtitles replace eform
esttab I Ip J Jp K L using "tableS4.csv", se(3) nogaps b(3) aic bic nonotes label star(* 0.10 ** 0.05 *** 0.01) nomtitles replace eform



********************************************************************************
// Models M, N, O, P, Q, R, and S (Table S5)
********************************************************************************

eststo M: stcox v2x_polyarchy2019 wbgi_gee2019 i.region $controls h_polcon5 if daterange, vce(robust) strata(date_firstcase)
eststo N: stcox v2x_polyarchy2019 wbgi_gee2019 i.region $controls dpi_presidentialism if daterange, vce(robust) strata(date_firstcase)
eststo O: stcox v2x_polyarchy2019 wbgi_gee2019 i.region $controls v2x_feduni if daterange, vce(robust) strata(date_firstcase)

stset time, failure(schoolsclosed_nat) id(ccode) scale(1)
eststo P:  stcox v2x_polyarchy2019 wbgi_gee2019 i.region $controls if daterange, vce(robust) strata(date_firstcase)
eststo Q:  stcox v2x_polyarchy2019 wbgi_gee2019 i.region $controls $controls2 if daterange, vce(robust) strata(date_firstcase)
stset time, failure(schoolsclosed) id(ccode) scale(1)

eststo R: stcox v2x_polyarchy2019 wbgi_gee2019 weekend i.region alt_ln_cases1 reg_any_case if daterange, vce(robust) strata(alt_date_firstcase)
eststo S: stcox v2x_polyarchy2019 wbgi_gee2019 weekend i.region alt_ln_cases1 reg_any_case $controls2 if daterange, vce(robust) strata(alt_date_firstcase)

esttab M N O P Q R S, se(3) nogaps b(3) aic bic nonotes label star(* 0.10 ** 0.05 *** 0.01) nomtitles replace eform
esttab M N O P Q R S using "tableS5.csv", se(3) nogaps b(3) aic bic nonotes label star(* 0.10 ** 0.05 *** 0.01) nomtitles replace eform



********************************************************************************
// Additional robustness checks mentioned in footnotes
********************************************************************************

*** Check correlations mentioned in footnote 2
corr v2x_polyarchy2019 h_polcon5 
corr v2x_polyarchy2019 dpi_presidentialism 
corr v2x_polyarchy2019 v2x_feduni


*** Check robustness to other stratifications mentioned in footnote 5

* Stratification on week of first case
gen week_firstcase = wofd(date_firstcase)
stcox v2x_polyarchy2019 wbgi_gee2019 i.region $controls if daterange, vce(robust) strata(week_firstcase) 
stcox v2x_polyarchy2019 wbgi_gee2019 i.region $controls $controls2 if daterange, vce(robust) strata(week_firstcase) 

* No stratification
stcox v2x_polyarchy2019 wbgi_gee2019 i.region $controls if daterange, vce(robust) 
stcox v2x_polyarchy2019 wbgi_gee2019 i.region $controls $controls2 if daterange, vce(robust) 



********************************************************************************
// Analyses of other containment and health measures (Figure 6)
********************************************************************************

*** Other domestic containment measures

* Baseline version

label var wbgi_gee2019 `" "Government" "Effectiveness" "'
est clear
foreach v in c2_workplaceclosing c3_cancelpublicevents c4_restgatherings c5_closepublictransport c6_stayathome c7_restmovement{
gen rec_`v' = `v'
recode rec_`v' (0=0)(1/3=1)
stset time, failure(rec_`v') id(ccode) scale(1)
bysort ccode: egen md_`v' = min(date) if inlist(`v',1,2,3)
bysort ccode: egen datef_`v' = max(md_`v')
eststo `v': stcox v2x_polyarchy2019 wbgi_gee2019 i.region $controls  if `v' != . , vce(robust) strata(date_firstcase)
drop rec_`v' md_`v' datef_`v'
}

coefplot (c2_workplaceclosing, label(Close Workplaces) mcolor(black) ciopts(lcolor(black))) (c3_cancelpublicevents, label(Cancel Public Events) mcolor(black) ciopts(lcolor(black))) (c4_restgatherings, label(Limits on Gatherings) mcolor(black) ciopts(lcolor(black))) (c5_closepublictransport, label(Close Public Transport) mcolor(black) ciopts(lcolor(black))) (c6_stayathome, label(Stay-at-Home Policy) mcolor(black) ciopts(lcolor(black))) (c7_restmovement, label(Restricted Movement) mcolor(black) ciopts(lcolor(black))), legend(region(lstyle(none)))  title("A (Larger Sample)") eform(*) keep(v2x_polyarchy2019 wbgi_gee2019) xline(1) scheme(s1mono) xscale(range(0,2.5)) xlabel(0(0.5)2.5) scale(0.85) name(containmentA, replace)

* Version including additional controls

foreach v in c2_workplaceclosing c3_cancelpublicevents c4_restgatherings c5_closepublictransport c6_stayathome c7_restmovement{
gen rec_`v' = `v'
recode rec_`v' (0=0)(1/3=1)
stset time, failure(rec_`v') id(ccode) scale(1)
bysort ccode: egen md_`v' = min(date) if inlist(`v',1,2,3)
bysort ccode: egen datef_`v' = max(md_`v')
eststo `v': stcox v2x_polyarchy2019 wbgi_gee2019 i.region $controls $controls2 if `v' != . , vce(robust) strata(date_firstcase)
drop rec_`v' md_`v' datef_`v'
}

coefplot (c2_workplaceclosing, label(Close Workplaces) mcolor(black) ciopts(lcolor(black))) (c3_cancelpublicevents, label(Cancel Public Events) mcolor(black) ciopts(lcolor(black))) (c4_restgatherings, label(Limits on Gatherings) mcolor(black) ciopts(lcolor(black))) (c5_closepublictransport, label(Close Public Transport) mcolor(black) ciopts(lcolor(black))) (c6_stayathome, label(Stay-at-Home Policy) mcolor(black) ciopts(lcolor(black))) (c7_restmovement, label(Restricted Movement) mcolor(black) ciopts(lcolor(black))), legend(region(lstyle(none)))  title("B (Additional Controls)") eform(*) keep(v2x_polyarchy2019 wbgi_gee2019) xline(1) scheme(s1mono) xscale(range(0,2.5)) xlabel(0(0.5)2.5) scale(0.85) name(containmentB, replace)

*** Less-disruptive health measures

* Baseline version

foreach v in h1_publicinfo h2_testingpolicy h3_contacttracing {
gen req_`v' = `v'
recode req_`v' (0/1=0)(2/3=1)
stset time, failure(req_`v') id(ccode) scale(1)
bysort ccode: egen md_`v' = min(date) if inlist(`v',1,2,3)
bysort ccode: egen datef_`v' = max(md_`v')
eststo `v': stcox v2x_polyarchy2019 wbgi_gee2019 i.region $controls if `v' != . , vce(robust) strata(date_firstcase)
drop req_`v' md_`v' datef_`v'
}

foreach v in h8_protectionofelderly{
gen rec_`v' = `v'
recode rec_`v' (0=0)(1/3=1)
stset time, failure(rec_`v') id(ccode) scale(1)
bysort ccode: egen md_`v' = min(date) if inlist(`v',1,2,3)
bysort ccode: egen datef_`v' = max(md_`v')
eststo `v': stcox v2x_polyarchy2019 wbgi_gee2019 i.region $controls if `v' != . , vce(robust) strata(date_firstcase)
drop rec_`v' md_`v' datef_`v'
}

coefplot (h1_publicinfo, label(Public Info Campaign) mcolor(gray) ciopts(lcolor(gray))) (h2_testingpolicy, label(Comprehensive Testing) mcolor(gray) ciopts(lcolor(gray))) (h3_contacttracing, label(Comprehensive Tracing) mcolor(gray) ciopts(lcolor(gray))) (h8_protectionofelderly, label(Isolation of Elderly) mcolor(gray) ciopts(lcolor(gray))), legend(region(lstyle(none)))  title("A (Larger Sample)") eform(*) keep(v2x_polyarchy2019 wbgi_gee2019) xline(1) xscale(range(0,2.5)) xlabel(0(0.5)2.5) scheme(s1mono) name(healthA, replace) scale(0.85)

* Version including additional controls

foreach v in h1_publicinfo h2_testingpolicy h3_contacttracing {
gen req_`v' = `v'
recode req_`v' (0/1=0)(2/3=1)
stset time, failure(req_`v') id(ccode) scale(1)
bysort ccode: egen md_`v' = min(date) if inlist(`v',1,2,3)
bysort ccode: egen datef_`v' = max(md_`v')
eststo `v': stcox v2x_polyarchy2019 wbgi_gee2019 i.region $controls $controls2 if `v' != . , vce(robust) strata(date_firstcase)
drop req_`v' md_`v' datef_`v'
}

foreach v in h8_protectionofelderly{
gen rec_`v' = `v'
recode rec_`v' (0=0)(1/3=1)
stset time, failure(rec_`v') id(ccode) scale(1)
bysort ccode: egen md_`v' = min(date) if inlist(`v',1,2,3)
bysort ccode: egen datef_`v' = max(md_`v')
eststo `v': stcox v2x_polyarchy2019 wbgi_gee2019 i.region $controls $controls2 if `v' != . , vce(robust) strata(date_firstcase)
drop rec_`v' md_`v' datef_`v'
}

coefplot (h1_publicinfo, label(Public Info Campaign) mcolor(gray) ciopts(lcolor(gray))) (h2_testingpolicy, label(Comprehensive Testing) mcolor(gray) ciopts(lcolor(gray))) (h3_contacttracing, label(Comprehensive Tracing) mcolor(gray) ciopts(lcolor(gray))) (h8_protectionofelderly, label(Isolation of Elderly) mcolor(gray) ciopts(lcolor(gray))), legend(region(lstyle(none))) title("B (Additional Controls)") eform(*) keep(v2x_polyarchy2019 wbgi_gee2019) xline(1) xscale(range(0,2.5)) xlabel(0(0.5)2.5) scheme(s1mono) name(healthB, replace) scale(0.85) 


graph combine containmentA containmentB, xsize(6) ysize(3) name(containment, replace) title("Other Domestic Containment Measures", size(medium)) 
graph combine healthA healthB, xsize(6) ysize(3) name(health, replace) title("Less-Disruptive Health Measures" , size(medium))
graph combine containment health, col(1) xsize(6) ysize(6) scale(0.87)
graph export "fig6.pdf", as(pdf) replace

* Set schoolsclosed as failure event again 
stset time, failure(schoolsclosed) id(ccode) scale(1)



********************************************************************************
// 170 versions of Models A and B, excluding one country at a time (Figure S4)
********************************************************************************

preserve

matrix coefsJK = J(173, 14, .)
matrix coln coefsJK = A_D_b A_D_ll95 A_D_ul95 A_SC_b A_SC_ll95 A_SC_ul95 A_excl B_D_b B_D_ll95 B_D_ul95 B_SC_b B_SC_ll95 B_SC_ul95 B_excl
local i=1
levelsof ccode, clean local(ccodes)
foreach excode of local ccodes {
quietly: stcox  v2x_polyarchy2019 wbgi_gee2019 i.region $controls if ccode != `excode' , vce(robust) strata(date_firstcase) 
matrix b = r(table)
matrix coefsJK[`i',1] = b[1,1]
matrix coefsJK[`i',2] = b[5,1]
matrix coefsJK[`i',3] = b[6,1]
matrix coefsJK[`i',4] = b[1,2]
matrix coefsJK[`i',5] = b[5,2]
matrix coefsJK[`i',6] = b[6,2]
matrix coefsJK[`i',7] = `excode'
quietly: stcox v2x_polyarchy2019 wbgi_gee2019 i.region $controls $controls2 if ccode != `excode' , vce(robust) strata(date_firstcase) 
matrix b = r(table)
matrix coefsJK[`i',8] = b[1,1]
matrix coefsJK[`i',9] = b[5,1]
matrix coefsJK[`i',10] = b[6,1]
matrix coefsJK[`i',11] = b[1,2]
matrix coefsJK[`i',12] = b[5,2]
matrix coefsJK[`i',13] = b[6,2]
matrix coefsJK[`i',14] = `excode'
local i=`i'+1
}

svmat2 coefsJK, rname(varname) name(col)
keep A_D_b A_D_ll95 A_D_ul95 A_SC_b A_SC_ll95 A_SC_ul95 A_excl B_D_b B_D_ll95 B_D_ul95 B_SC_b B_SC_ll95 B_SC_ul95 B_excl
keep if _n <=173

* Omit China, Hong Kong and Mongolia
drop if A_excl == 156
drop if A_excl == 344
drop if A_excl == 496	
gen n = _n
label var n "Exclusion number (of 170)"

label var A_D_b "Democracy"
label var A_D_ll95 "Democracy (95% CI)"
label var A_D_ul95 "Democracy (95% CI)"
label var A_SC_b "Goverment Effectiveness"
label var A_SC_ll95 "Govt Effectiveness (95% CI)"
label var A_SC_ul95 "Govt Effectiveness (95% CI)"
line A_D_b A_D_ll95 A_D_ul95 A_SC_b A_SC_ll95 A_SC_ul95 n, yline(1, lpattern(dash)) lpattern(solid dot dot solid dot dot) lcolor(black black black gray gray gray ) xscale(range(0,170)) xlabel(0(34)170) title("A (Larger Sample)") legend(order(1 4) cols(1) region(lstyle(none))) name(a, replace) yscale(range(.25,1.75)) ylabel(.25(0.25)1.75)

label var B_D_b "Democracy"
label var B_D_ll95 "Democracy (95% CI)"
label var B_D_ul95 "Democracy (95% CI)"
label var B_SC_b "Government Effectiveness"
label var B_SC_ll95 "Govt Effectiveness (95% CI)"
label var B_SC_ul95 "Govt Effectiveness (95% CI)"
line B_D_b B_D_ll95 B_D_ul95 B_SC_b B_SC_ll95 B_SC_ul95 n, yline(1, lpattern(dash)) lpattern(solid dot dot solid dot dot)  lcolor(black black black gray gray gray) xscale(range(0,170)) xlabel(0(34)170) title("B (Additional Controls)") legend(order(1 4 ) cols(1) region(lstyle(none))) name(b, replace) yscale(range(.25,1.75)) ylabel(.25(0.25)1.75)

graph combine a b, iscale(0.75) l1title("Hazard ratio estimates and 95% CIs")
graph export "figS4.pdf", as(pdf) replace

restore



********************************************************************************
// Country-wise list of key variables (Table S6)
********************************************************************************

format date_firstcase date_firstclosure date_firstclosure_nat %tdMon_dd
clonevar print_date_firstclosure = date_firstclosure
replace print_date_firstclosure = mdy(4,7,2020) if country == "Nicaragua"
* Copy data manually from data browser:
br country date_firstcase date_firstclosure date_firstclosure_nat cases v2x_polyarchy2019 wbgi_gee2019 if date == print_date_firstclosure
drop print_date_firstclosure



********************************************************************************
// Export eight datasets for Cox ED analysis in R (see separate R markdown file)
********************************************************************************

preserve
drop if date < mdy(01,28,2020)
drop if date > mdy(04,07,2020)
format date_firstcase %tg
gen date2 = date
format date2 %tg
gen date1 = date2-1

* A+B (V-Dem + WGI)
xtreg schoolsclosed v2x_polyarchy2019 wbgi_gee2019  $controls, re 
sum wbgi_gee2019 if e(sample), de
gen lowSC = r(p25)
gen highSC = r(p75)
sum v2x_polyarchy2019 if e(sample), de
gen lowD = r(p25)
gen highD = r(p75)
stcox  v2x_polyarchy2019 wbgi_gee2019  $controls, 
export delimited schoolsclosed time ccode date date1 date2 date_firstcase v2x_polyarchy2019 fh_ipolity2 weekend wbgi_gee2019 ln_cases1 _st region $controls $controls2 low* high* using "JPP_covid19_2022_CoxED_A.csv" if e(sample), replace
stcox  v2x_polyarchy2019 wbgi_gee2019  $controls $controls2, 
export delimited schoolsclosed time ccode date date1 date2 date_firstcase v2x_polyarchy2019 fh_ipolity2 weekend wbgi_gee2019 ln_cases1 _st region $controls $controls2 low* high* using "JPP_covid19_2022_CoxED_B.csv" if e(sample), replace

*C+D (Vanhanen + WGI)
xtreg schoolsclosed van_comp van_part v2x_polyarchy2019 wbgi_gee2019  $controls, re 
sum van_part if e(sample), de
gen lowP = r(p25)
gen highP = r(p75)
sum van_comp if e(sample), de
gen lowC = r(p25)
gen highC = r(p75)
eststo: stcox  van_comp van_part wbgi_gee2019 i.region $controls , vce(robust) strata(date_firstcase)
export delimited schoolsclosed van_comp van_part time ccode date date1 date2 date_firstcase v2x_polyarchy2019 fh_ipolity2 weekend wbgi_gee2019 ln_cases1 _st region $controls  low* high* using "JPP_covid19_2022_CoxED_C.csv" if e(sample), replace
eststo: stcox  van_comp van_part wbgi_gee2019 i.region $controls $controls2 , vce(robust) strata(date_firstcase)
export delimited schoolsclosed van_comp van_part time ccode date date1 date2 date_firstcase v2x_polyarchy2019 fh_ipolity2 weekend wbgi_gee2019 ln_cases1 _st region $controls $controls2 low* high* using "JPP_covid19_2022_CoxED_D.csv" if e(sample), replace

*E+F (FH/Polity + ICRG)
xtreg schoolsclosed van_comp van_part v2x_polyarchy2019 wbgi_gee2019  $controls, re 
sum fh_ipolity2  if e(sample), de
gen lowFPDem = r(p25)
gen highFPDem = r(p75)
sum icrg_qog if e(sample), de
gen lowQOG = r(p25)
gen highQOG = r(p75)
eststo: stcox  fh_ipolity2 icrg_qog i.region $controls , vce(robust) strata(date_firstcase)
export delimited schoolsclosed fh_ipolity2 icrg_qog time ccode date date1 date2 date_firstcase v2x_polyarchy2019 fh_ipolity2 weekend wbgi_gee2019 ln_cases1 _st region $controls  low* high* using "JPP_covid19_2022_CoxED_E.csv" if e(sample), replace
eststo: stcox  fh_ipolity2 icrg_qog  i.region $controls $controls2 , vce(robust) strata(date_firstcase)
export delimited schoolsclosed fh_ipolity2 icrg_qog time ccode date date1 date2 date_firstcase v2x_polyarchy2019 fh_ipolity2 weekend wbgi_gee2019 ln_cases1 _st region $controls $controls2 low* high* using "JPP_covid19_2022_CoxED_F.csv" if e(sample), replace

*G+H (Polity PARCOMP/PARREG + ICRG)

xtreg schoolsclosed parreg parcomp icrg_qog v2x_polyarchy2019 wbgi_gee2019 $controls, re 
sum parreg if e(sample), de
gen lowPolP = r(p25)
gen highPolP = r(p75)
sum parcomp if e(sample), de
gen lowPolC = r(p25)
gen highPolC = r(p75)
eststo: stcox  parreg parcomp icrg_qog i.region $controls , vce(robust) strata(date_firstcase)
export delimited schoolsclosed parreg parcomp icrg_qog time ccode date date1 date2 date_firstcase v2x_polyarchy2019 fh_ipolity2 weekend wbgi_gee2019 ln_cases1 _st region $controls  low* high* using "JPP_covid19_2022_CoxED_G.csv" if e(sample), replace
eststo: stcox  parreg parcomp icrg_qog i.region $controls $controls2 , vce(robust) strata(date_firstcase)
export delimited schoolsclosed parreg parcomp icrg_qog time ccode date date1 date2 date_firstcase v2x_polyarchy2019 fh_ipolity2 weekend wbgi_gee2019 ln_cases1 _st region $controls $controls2 low* high* using "JPP_covid19_2022_CoxED_H.csv" if e(sample), replace

restore



********************************************************************************
// Manually enter and plot Cox ED output from R (Figures 4, 5, S2, and S3)
********************************************************************************

*** Coefplots based on results from R. Minor differences in the
*** confidence intervals can occur due to the bootstrapping procedure

*** Figure 4 

matrix coefsA = J(6, 3, .)
matrix coln coefsA = effect ll95 ul95
matrix rown coefsA = dem govt demlowgovt demhighgovt govtlowdem govthighdem
matrix coefsA[1,1] = -10.999
matrix coefsA[1,2] = -14.373
matrix coefsA[1,3] = -7.625
matrix coefsA[2,1] = 7.858  
matrix coefsA[2,2] = 5.459
matrix coefsA[2,3] = 10.258
matrix coefsA[3,1] = -16.283
matrix coefsA[3,2] = -21.564
matrix coefsA[3,3] = -11.002
matrix coefsA[4,1] = -6.675
matrix coefsA[4,2] = -8.997
matrix coefsA[4,3] = -4.354
matrix coefsA[5,1] = 4.329
matrix coefsA[5,2] = 2.815
matrix coefsA[5,3] = 5.843
matrix coefsA[6,1] = 13.936
matrix coefsA[6,2] = 9.492
matrix coefsA[6,3] = 18.381
matrix list coefsA

matrix coefsB = J(6, 3, .)
matrix coln coefsB = effect ll95 ul95
matrix rown coefsB = dem govt demlowgovt demhighgovt govtlowdem govthighdem
matrix coefsB[1,1] = -8.220
matrix coefsB[1,2] = -9.762
matrix coefsB[1,3] = -6.679
matrix coefsB[2,1] = 8.275
matrix coefsB[2,2] = 6.613
matrix coefsB[2,3] = 9.937
matrix coefsB[3,1] = -11.852
matrix coefsB[3,2] = -14.064
matrix coefsB[3,3] = -9.641
matrix coefsB[4,1] = -5.152
matrix coefsB[4,2] = -6.342
matrix coefsB[4,3] = -3.961
matrix coefsB[5,1] = 6.006
matrix coefsB[5,2] = 4.654
matrix coefsB[5,3] = 7.357
matrix coefsB[6,1] = 12.706
matrix coefsB[6,2] = 10.245
matrix coefsB[6,3] = 15.168
matrix list coefsB

mylabels -25(5)25, myscale(@/1) local(myla)
coefplot (matrix(coefsA[,1]), ci((coefsA[,2] coefsA[,3])) label(A & A + (Interaction)) msymbol(circle_hollow) mcolor(black) ciopts(lcolor(black))) (matrix(coefsB[,1]), ci((coefsB[,2] coefsB[,3])) label(B & B + (Interaction)) msymbol(circle) mcolor(black) ciopts(lcolor(black))) , headings(dem = "{bf:Baseline Models}" demlowgovt = "{bf:Interaction Models}") coeflabels(dem = "Democracy" govt = "Government Effectiveness" demlowgovt = `""Democracy at" "Low Government Effectiveness""' demhighgovt  = `""Democracy at" "High Government Effectiveness""' govtlowdem  = `""Government Effectiveness at" "Low Level of Democracy""' govthighdem = `""Government Effectiveness at" "High Level of Democracy""') xline(0)  xlab(`myla') xtitle("Difference in Expected Duration Between" "Observations at 75th and 25th Percentile") legend(region(lstyle(none)))  xsize(10) ysize(6) grid(between) scale(0.9) aspectratio(0.8)
graph export "fig4.pdf", as(pdf) replace


*** Figure 5 

matrix coefsC = J(3, 3, .)
matrix coln coefsC = effect ll95 ul95
matrix rown coefsC = "Closed Authoritarian Regime" "Electoral Authoritarian Regime" "Democratic Regime"
matrix coefsC[1,1] = 0
matrix coefsC[1,2] = 0
matrix coefsC[1,3] = 0
matrix coefsC[2,1] = -1.263
matrix coefsC[2,2] = -1.636
matrix coefsC[2,3] = -0.891
matrix coefsC[3,1] = -7.093  
matrix coefsC[3,2] = -9.236
matrix coefsC[3,3] = -4.949
matrix list coefsC

matrix coefsD = J(3, 3, .)
matrix coln coefsD = effect ll95 ul95
matrix rown coefsD = "Closed Authoritarian Regime" "Electoral Authoritarian Regime" "Democratic Regime"
matrix coefsD[1,1] = 0
matrix coefsD[1,2] = 0
matrix coefsD[1,3] = 0
matrix coefsD[2,1] = 0.854
matrix coefsD[2,2] = 0.702
matrix coefsD[2,3] = 1.005
matrix coefsD[3,1] = -3.705
matrix coefsD[3,2] = -4.322
matrix coefsD[3,3] = -3.087
matrix list coefsD

mylabels -10(2)10, myscale(@/1) local(myla)
coefplot (matrix(coefsC[,1]), ci((coefsC[,2] coefsC[,3])) label(C (Larger Sample)) msymbol(diamond_hollow) mcolor(gray) ciopts(lcolor(gray))) (matrix(coefsD[,1]), ci((coefsD[,2] coefsD[,3])) label(D (Additional Controls)) msymbol(diamond) mcolor(gray) ciopts(lcolor(gray))), xline(0) xlab(`myla') xtitle("Difference in Expected Duration Compared" "to a Closed Authoritarian Regime") legend(region(lstyle(none)))  xsize(10) ysize(6)  grid(between) scale(0.9) aspectratio(0.8)
graph export "fig5.pdf", as(pdf) replace


*** Figure S2

matrix coefsE = J(6, 3, .)
matrix coln coefsE = effect ll95 ul95
matrix rown coefsE = dem govt demlowgovt demhighgovt govtlowdem govthighdem
matrix coefsE[1,1] =  -14.121
matrix coefsE[1,2] = -21.171
matrix coefsE[1,3] = -7.071
matrix coefsE[2,1] = 8.283  
matrix coefsE[2,2] = 4.144
matrix coefsE[2,3] = 12.422
matrix coefsE[3,1] = -16.896 
matrix coefsE[3,2] = -25.339
matrix coefsE[3,3] = -8.453
matrix coefsE[4,1] = -9.961
matrix coefsE[4,2] = -15.114
matrix coefsE[4,3] = -4.808
matrix coefsE[5,1] = 3.579
matrix coefsE[5,2] = 1.628 
matrix coefsE[5,3] = 5.531
matrix coefsE[6,1] = 10.514 
matrix coefsE[6,2] = 4.940
matrix coefsE[6,3] = 16.089
matrix list coefsE

matrix coefsF = J(6, 3, .)
matrix coln coefsF = effect ll95 ul95
matrix rown coefsF = dem govt demlowgovt demhighgovt govtlowdem govthighdem
matrix coefsF[1,1] = -9.432
matrix coefsF[1,2] = -14.045
matrix coefsF[1,3] = -4.819
matrix coefsF[2,1] = 5.993 
matrix coefsF[2,2] = 2.675
matrix coefsF[2,3] = 9.311
matrix coefsF[3,1] = -11.474
matrix coefsF[3,2] = -17.143
matrix coefsF[3,3] = -5.805
matrix coefsF[4,1] = -7.284
matrix coefsF[4,2] = -10.970
matrix coefsF[4,3] = -3.598
matrix coefsF[5,1] = 3.141
matrix coefsF[5,2] = 1.207
matrix coefsF[5,3] = 5.075
matrix coefsF[6,1] = 7.331
matrix coefsF[6,2] = 3.502
matrix coefsF[6,3] = 11.160
matrix list coefsF

mylabels -25(5)25, myscale(@/1) local(myla)
coefplot (matrix(coefsE[,1]), ci((coefsE[,2] coefsE[,3])) label(E & E + (Interaction)) msymbol(circle_hollow) mcolor(black) ciopts(lcolor(black))) (matrix(coefsF[,1]), ci((coefsF[,2] coefsF[,3])) label(F & F + (Interaction)) msymbol(circle) mcolor(black) ciopts(lcolor(black))), headings(dem = "{bf:Baseline Models}" demlowgovt = "{bf:Interaction Models}") coeflabels(dem = "Democracy (FH/Polity)" govt = "Quality of Government (ICRG)" demlowgovt = `""Democracy at" "Low Quality of Government""' demhighgovt  = `""Democracy at" "High Quality of Government""' govtlowdem  = `""Quality of Government at" "Low Level of Democracy""' govthighdem = `""Quality of Government at" "High Level of Democracy""') xline(0) xlab(`myla') xtitle("Difference in Expected Duration Between" "Observations at 75th and 25th Percentile") legend(region(lstyle(none)))  xsize(10) ysize(6) grid(between) scale(0.9) aspectratio(0.8)
graph export "figS2.pdf", as(pdf) replace


*** Figure S3

matrix coefsG = J(3, 3, .)
matrix coln coefsG = effect ll95 ul95
matrix rown coefsG = "Closed Authoritarian Regime" "Electoral Authoritarian Regime" "Democratic Regime"
matrix coefsG[1,1] = 0
matrix coefsG[1,2] = 0
matrix coefsG[1,3] = 0
matrix coefsG[2,1] = 1.045
matrix coefsG[2,2] = 0.401
matrix coefsG[2,3] = 1.688
matrix coefsG[3,1] = -13.304 
matrix coefsG[3,2] = -20.733
matrix coefsG[3,3] = -5.874
matrix list coefsG

matrix coefsH = J(3, 3, .)
matrix coln coefsH = effect ll95 ul95
matrix rown coefsH = "Closed Authoritarian Regime" "Electoral Authoritarian Regime" "Democratic Regime"
matrix coefsH[1,1] = 0
matrix coefsH[1,2] = 0
matrix coefsH[1,3] = 0
matrix coefsH[2,1] = 0.608
matrix coefsH[2,2] = 0.161
matrix coefsH[2,3] = 1.054
matrix coefsH[3,1] = -9.854
matrix coefsH[3,2] = -15.827
matrix coefsH[3,3] = -3.882
matrix list coefsH


mylabels -25(5)25, myscale(@/1) local(myla)
coefplot (matrix(coefsG[,1]), ci((coefsG[,2] coefsG[,3])) label(G (Larger Sample)) msymbol(diamond_hollow) mcolor(gray) ciopts(lcolor(gray))) (matrix(coefsH[,1]), ci((coefsH[,2] coefsH[,3])) label(H (Additional Controls)) msymbol(diamond) mcolor(gray) ciopts(lcolor(gray))), xline(0) xlab(`myla') xtitle("Difference in Expected Duration Compared" "to a Closed Authoritarian Regime") legend(region(lstyle(none)))  xsize(10) ysize(6)  grid(between) scale(0.9) aspectratio(0.8) 
graph export "figS3.pdf", as(pdf) replace



********************************************************************************
// End
********************************************************************************
